﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;
using LatoolNet;

namespace NumericalSolution {
  class Diffusion2DADI {

    private double m_deltax;
    private double m_deltay;
    private double m_deltat;
    private double m_alpha;
    private int m_nx;
    private int m_ny;
    private double m_boundaryCondition = 0.0;

    public int NumberOfNodesInX {
      get { return m_nx; }
      set { m_nx = value; }
    }

    public int NumberOfNodesInY {
      get { return m_ny; }
      set { m_ny = value; }
    }

    public double DeltaX {
      get { return m_deltax; }
      set { m_deltax = value; }
    }

    public double DeltaY {
      get { return m_deltay; }
      set { m_deltay = value; }
    }

    public double DeltaT {
      get { return m_deltat; }
      set { m_deltat = value; }
    }

    public double Alpha {
      get { return m_alpha; }
      set { m_alpha = value; }
    }

    public void NextStep(double[,] u) {

      double[,] u_half = new double[m_nx, m_ny];

      for (int j = 1; j < m_ny - 1; j++) {
        StepI(u, u_half, j);
      }

      for (int i = 1; i < m_nx - 1; i++) {
        StepII(u, u_half, i);
      }

    }

    private void StepI(double[,] u, double[,] u_half, int j) {

      Matrix mat = new Matrix(m_nx, m_nx, 3);
      Matrix vec = new Matrix(m_nx, 1);

      mat[0, 0] = 1;
      vec[0, 0] = m_boundaryCondition;

      double rx;
      double ry;

      for (int i = 1; i < m_nx - 1; i++) {

        rx = (m_alpha * m_deltat) / (m_deltax * m_deltax);
        ry = (m_alpha * m_deltat) / (m_deltay * m_deltay);

        mat[i, i - 1] = -rx;
        mat[i, i] = (1 + 2 * rx);
        mat[i, i + 1] = -rx;
        vec[i, 0] = ry * u[i, j - 1] 
                  + (1 - 2 * ry) * u[i, j] 
                  + ry * u[i, j + 1];
      }

      mat[m_nx - 1, m_nx - 1] = 1;
      vec[m_nx - 1, 0] = m_boundaryCondition;

      //LUFactorization.Solve(mat, vec); 
      LinearEquation.Solve(mat, vec);

      for (int i = 0; i < m_nx; i++) {
        u_half[i, j] = vec[i, 0];
      }

      mat.Dispose();
      vec.Dispose();

    }

    private void StepII(double[,] u, double[,] u_half, int i) {

      Matrix mat = new Matrix(m_ny, m_ny, 3);
      Matrix vec = new Matrix(m_ny, 1);

      mat[0, 0] = 1;
      vec[0, 0] = m_boundaryCondition;

      double rx;
      double ry;

      for (int j = 1; j < m_ny - 1; j++) {

        rx = (m_alpha * m_deltat) / (m_deltax * m_deltax);
        ry = (m_alpha * m_deltat) / (m_deltay * m_deltay);

        mat[j, j - 1] = -ry;
        mat[j, j] = (1 + 2 * ry);
        mat[j, j + 1] = -ry;
        vec[j, 0] = rx * u_half[i - 1, j] 
                  + (1 - 2 * rx) * u_half[i, j] 
                  + rx * u_half[i + 1, j];
      }

      mat[m_ny - 1, m_ny - 1] = 1;
      vec[m_ny - 1, 0] = m_boundaryCondition;

      //LUFactorization.Solve(mat, vec);
      LinearEquation.Solve(mat, vec);

      for (int j = 0; j < m_ny; j++) {
        u[i, j] = vec[j, 0];
      }

      mat.Dispose();
      vec.Dispose();
    }
  }
}
